!> author: 左志华
!> date: 2022-07-07
!>
!> Time integration <br>
!> 时间步迭代
!> @note experimental code
module sph_time_integration

    use sph_kinds, only: rk
    use sph_project_type, only: project_t
    use sph_region, only: region_t
    use sph_save_pif, only: save_PIF, nstep
    use sph_leap_frog, only: leap_frog_init, leap_frog_calc
    use sph_navier_stokes_solver, only: navier_stokes_solver, nsave
    use sph_terminal, only: term
    use sph_env, only: has_rigid, dim, nums

contains

    !> Time integration <br>
    !> 时间步迭代
    impure subroutine time_integration(project, region, nstart, timesteps, time_start, dt)
        type(project_t), intent(inout) :: project   !! project configuration <br>
                                                    !! 项目
        type(region_t), intent(inout) :: region     !! region of particle macroscopic properties <br>
                                                    !! 粒子宏观变量域
        integer, intent(in) :: nstart               !! start index of time integration <br>
                                                    !! 开始时间步
        integer, intent(in) :: timesteps            !! number of time steps <br>
                                                    !! 需要计算的时间步数
        real(rk), intent(in) :: time_start          !! start time of time integration <br>
                                                    !! 开始时间
        real(rk), intent(in) :: dt                  !! time step size <br>
                                                    !! 时间步长
        integer :: i, nstart_
        real(rk) :: t

        ! @todo 龙格库塔法

        t = time_start
        nstart_ = nstart

        if (project%itime_integration == 1) then
            if (.not. project%go_on) then   ! 如果是第一次运算，则初始化蛙跳法，否则需要处理蛙跳法续算 @todo
                call leap_frog_init(region, nstart_, t, dt, navier_stokes_solver)
                nstart_ = nstart_ + 1
                t = t + dt
                if (mod(nstart_, nsave) == 0) &
                    call save_PIF(project%fid, region, dim, nums(4), t, nstep)
                project%go_on = .true.
                
                ! @todo rigid
                ! if (has_rigid) then
                !     call 
                
            end if
        end if

        do i = nstart_, nstart - 1 + timesteps
            t = t + dt
            call project%time_integration(region, i, t, dt, navier_stokes_solver)     ! 进行时间步迭代，并通过回调函数求解加速度
            if (mod(i, nsave) == 0) then
                call save_PIF(project%fid, region, dim, nums(4), t, nstep)
                call term%progress_bar('Time integration:', real(i)/real(timesteps))
            end if
        end do
        call term%progress_bar('Time integration:', 1.0)

    end subroutine time_integration

end module sph_time_integration
